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(3JT)i By numerical simulations of the 3d Ising spin glass we find evidence that spontaneous replica 

^ ' symmetry breaking theory and not the droplet model describes with good accuracy the equilibrium 

, behavior of the system. 

O 1 PACS numbers: 75.50.Lk, 05.50.+q, 64.60. Cn, 02.50.Ng, 02.60.Cb 



The behavior of disordered spin models at equilibrium 
is well understood in the framework of the mean field ap- 
proximation QJ^] . The main prediction of the mean field 
approach is the existence of a low temperature glassy 
phase. Such a phase is characterized by the existence 
of many different equilibrium states (spontaneous replica 
symmetry breaking, SRSB). On the other hand it is pos- 
sible to define a different consistent theory || by starting 
from the Migdal-Kadanoff approach. By using a common 
terminology we will refer in the following to this approach 
as to the droplet model. Here one expects the equilib- 
rium state to be unique (apart from global inversions in 
zero magnetic field) and that the most relevant excita- 
tions are obtained by reversing large domains of spins 
q ■ (the droplets). 

We have two different starting points. One is the infi- 
nite range approximation which leads to the replica sym- 
metry breaking picture and the other is the Kadanoff- 
Migdal approximation which leads to the droplet model. 
Although each of the two pictures is correct in its range 
of validity we have to establish which of the two qualita- 
tively describes the physics of the real three dimensional 
spin glasses. 

The main result of this work (which continues the in- 
vestigation started in Q, and follows a long series of 
Monte Carlo simulations of spin glass systems [pi) has 
been to gather new and strong evidence that in three di- 
mensions the SRSB picture (and not the droplet model) 
describes correctly what is observed in numerical simu- 
lations. 

Let us start by summarizing the evidence we will 
present in this note and the scheme of our reasoning (for 
a more detailed exposition of these and more data see ref. 
H ) . We will start by showing that the probability distri- 
bution of the overlap among two systems at equilibrium, 
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P(q), has a non-trivial structure. P(q ~ 0) is different 
from zero, and its shape does not depend on the volume 
size. We will analyze (following a suggestion contained in 
the third reference of 0) sample to sample fluctuations 
of the spin-glass susceptibility, and find out that they 
are incompatible with the droplet model, while their size 
is very well explained (even in a quantitative manner) 
by SRSB theory. In order to show that the structure of 
the different equilibrium states is not compatible with a 
droplet structure we will compute and analyze equal time 
correlation functions. From this analysis we deduce the 
existence of many equilibrium states that cannot be de- 
scribed by a droplet like structure. On the contrary we 
will show that even at a quantitative level SRSB theory 
explains very well the numerical data. 

Further evidence about the inadequacy of the droplet 
model to describe the 3d spin glasses and support to a 
SRSB mechanism will be provided by analyzing the dis- 
tribution of overlaps of boxes of side R, and by discussing 
the behavior of the box overlap Binder parameter. 

The model we will consider most is defined by the sim- 
ple Edwards- Anderson Hamiltonian on a 3d simple cubic 
lattice H = — k} a iJiM a ji where the sum runs over 
nearest neighbor couples of sites. The quenched disor- 
dered couplings J are distributed according to a Gaus- 
sian law. A study of the overlap susceptibility and of 
the Binder cumulants shows that (under the a-priori as- 
sumption the existence of a phase transition at a non-zero 
temperature with a power law divergence) the transition 
is located at T « 1. In order to check universality of 
our results we have also studied a model || with integer 
J = ±1 variables, where each spin is coupled with equal 
strength to 26 neighboring sites (all the ones contained 
in a cube of 3 3 sites). The results we discuss here are 
confirmed by our findings about this second model. 

We have used an isotropic lattice of linear size L, and 
we have computed the probability distribution Pj(q) of 
the overlap q = V^ x ^2i i OiTi among two thermalized con- 
figurations cr and t in a box of volume V — L 3 . We have 
studied the behavior of the function P(q) averaged over 
a large number of realizations of the quenched disordered 
couplings J (i.e. the average over the J random variables 
of Pj(q)). We have used a maximum of 2560 samples for 
the smallest lattice sizes and a minimum of 512 samples 
for the largest sizes. It was already known (see for exam- 
ple and references therein) that it is non-trivial and 
it has a shape quite similar to the one predicted in the 
mean field model. Mainly thanks to the use of large com- 
puter resources (we have mainly used the APE parallel 
computer H , which turns out to be very effective for this 
kind of problems ||: we flip about 2 -10 s spins per second 
on the tower version of the machine ) and of the temper- 
ing (an annealing-like improved Monte Carlo technique 
introduced in |B|) we have been able to study systems of 
larger size than before (up to 14 3 ), and to bring them to 
thermal equilibrium quite deep in the cold phase. In this 
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case we have equilibrated the system up to distance 14. 
We will see that this information is complemented by our 
dynamical study, where we work on time scales on which 
we can equilibrate the system on distances up to order 6. 
This gives a good control over the fractal geometry of the 
typical excitations and of their boundaries. This is what 
we need in order to distinguish between SRSB theory or 
Migdal-Kadanoff droplets. 

The first crucial comment is that the shape of the func- 
tion P(q) is in our statistical precision size independent. 
What we stress is that all the non-zero plateau at low q 
values, down to q — 0, turns out to be size- independent. 
For example at T = 0.7 the Binder cumulant of q is prac- 
tically independent from the lattice size and it is equal to 
0.85 ± .01. This means that the system has a non-trivial 
structure of equilibrium states with a continuous distri- 
bution of the allowed overlap values. By using our mea- 
surements of equal time correlation functions we will ar- 
gue in the following that such states cannot be described 
by the droplet approach, while they have all features pre- 
dicted by the SRSB approach. 

In this note we do not answer a very important ques- 
tion, i.e. if in the infinite volume limit a low-temperature 
phase characterized by the existence of a non zero order 
parameter qEA exists. On the lattice volumes we are able 
to investigate the high q peak of the P(q) is very slowly 
shifting toward lower q values, even if, as we already said, 
the shape of the P(q) does not change. The extrapolation 
to the infinite volume limit looks in this case a very del- 
icate issue, and many potential systematic errors (even 
in the definition of the finite volume qEA) are involved. 
Here we will not address in detail this point, and assume 
that we are working in conditions where the system is ef- 
fectively frozen to a phase with a non-zero value of qEA- 
A possible scenario [i| of a correlation length diverging 
exponentially for T — > or of a Kosterlitz-Thouless like 
transition would be compatible with this approach, since 
on our finite lattice we would be measuring properties of 
a frozen system. It is also important to note that this 
ambiguity only concerns the behavior of the high q peak 
of the P(q) (which could tend to q — on very large lat- 
tices), while on the contrary the P(q) for small q values 
is non-trivial and does not depend on the lattice size. 

The agreement with mean field theory becomes quan- 
titative if we study sample to sample fluctuations. Mean 
field theory tells us how much the function Pj(q) for a 
given realization of the quenched disorder differs from 
the average. For example if we consider (q k )j = 
J dq P.]{q) q k , we have that in mean field 

1- 



(q k )j(q m )j = ^(q k )j {q m )j + ^{q k+m )j , (l) 

where by the overline we indicate the average over the 
quenched noise. We have verified that in the low T region 
this equality is very well satisfied. For example for k = 2 
and m — 2 at T = 0.7 and L ranging from 4 to 10 the 
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ratio of the l.h.s. to the r.h.s. of ([!]) is equal to 1.0 with 
an error never larger than 0.1. 

Strictly speaking the non-triviality of the function P(q) 
is not in violent contradiction with the droplet model. In 
the framework of the droplet approach it is always pos- 
sible to suppose that states where domains that take a 
finite part of the whole system are reversed have a finite 
probability. This hypothesis is however rather unnatural 
and it is definitely wrong in the Kadanoff-Migdal approx- 
imation. Moreover we have already seen that the ability 
of the SRSB theory to predict quantitatively the fluctu- 
ations of the function P(q) is remarkable. 

We will now falsify the possibility discussed in the last 
paragraph by considering the q — q correlation functions 
restricted to those pairs of configuration which have a 
small value of q. The analysis of such correlation func- 
tions, together with the non-triviality of the P(q), will 
constitute an ultimate test of the failure of the droplet 
model. 

More precisely we consider a system of side L and we 
define the relevant correlation function as 

C(x, L) ^v-^Y. °i+*n+x Pin) , (2) 

i 

where the brackets indicate the thermal average. The 
droplet model predicts that C(x, oo) goes to the constant 
value q%A f° r l ar § e x - I n the SRSB approach C(x, oo) oc 
|a;|~ A , where A in an appropriate exponent which has 
been computed in less than 6 dimensions for the q = 
correlation functions [|o| . 

We have studied this problem by considering large sys- 
tems, with L = 64. We have ran numerical simulations 
starting from two random configurations selected inde- 
pendently (for 4 realizations of the quenched couplings). 
We have verified that q 2 stays small in the whole run so 
that the difference in the initial configurations, for not 
too large times, enforces the condition q w 0. Eventually 
in a finite system global equilibrium will be reached and 
q will become of order 1 . However if we let L — ► oo first 
we can use this approach to study the equilibrium value 
of the correlation function with the constraint of having 
zero overlap. 

In order to do that we consider the time dependent 
equal time correlation function at time t 

G(x,t) = V^'^^i+xTi+x <TiTi)t , (3) 
i 

where the average is done at time t, i.e. after t Monte 
Carlo cycles after the random start. We find that for 
large times t the correlation function G(x, t) is essen- 
tially different from zero for distances not too larger than 
a dynamic correlation lengthy £(t) which increases (and 
maybe diverges) with time. Our numerical data are well 
represented with the functional form 
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G(x,t) = 




where we have defined £(T, t) = B(T) t x{T ^ . In the whole 
range of distances 1 < x < 8 for Monte Carlo times which 
range from 10 2 to 10 6 full lattice sweeps and a large range 
of temperatures T < T c (we have done measurements 
at different temperatures, down to T m i n ~ .3T C ) we get 
good fits. The exponents a and 5 are weakly dependent 
on T. For example at T = 0.70 we get the best values a = 
0.50 ± 0.02 and <5 = 1.48 ± 0.02. The correlation length 
exponent A(T) is approximately given by 0.16T. Such 
power law growth of the correlation length was already 
observed by Rieger fuJI . In order to study the limit t — » 
oo in a safe way it is even better to avoid global fits and 
to fit the data at fixed distance x as 

G(x,t) = G(x,oo) exp{-A(x) t~ x( - T A . (5) 

In this way the extrapolation to infinite time (with the 
self-implemented constraint of q — always satisfied) is 
performed in a very safe way. We plot in Fig. (Q) the 
correlations G(x, oo) (computed at T — 0.7) as a function 
of the distance in double logarithmic scale. 

We have also computed the same quantities by using 
a different temperature schedule. In this second numer- 
ical experiment we slowly cool down the system from 
T = 1.5 > T c to the final temperature. To perform the 
cooling we use a number of steps proportional to t, the 
waiting time we want to look the correlation function at. 
After that the system evolves at the fixed temperature 
of interest T for t more time steps before measurement. 
In this way one can obtain a much better equilibration. 
As matter of principle in this case one does not expect 
a pure power law but a combination of different powers 
generated by different temperature contributions. How- 
ever a fit similar to the previous one (|J) works very well 
with a slightly large value of of A. At T = 0.7 one obtains 
the results shown in Fig. (]]]). The data obtained with 
the two techniques behave in a very similar way. The 
t = oo data are well described by a power decay x~ a 
with a — 0.50 ± .03, as predicted by the replica theory 
and in variance with the droplet model predictions. 

The value of the correlation function at distance 1 is 
particularly interesting. Indeed in the model with Gaus- 
sian quenched disorder one can easily prove by a simple 
integration by part that E = —/3(1 — C(l)), where E is 
the energy per link and C(x) is the overlap correlation 
function of the fully equilibrated system (i.e. summed 
over different ergodic components). Only if replica sym- 
metry is broken it differs by our correlation function G 
where the two replicas have been kept (by their own will) 
in two different ergodic components with zero mutual 
overlap. The energy can be computed with high accuracy 
and from its value we can deduce C(l). The equality 

C(1)=G(1) (6) 
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should be violated as soon as replica symmetry is broken, 
in the same way in which the equality E = — ^(1 — Q^a) 
is violated in the Sherrington-Kirkpatrick model at low 
temperature. 

The value of the energy per link is very well fitted by 
the form + At~ A( - T \ The exponent A(T) turns out 
to be quite large, i.e. we find A(T) ~ 0.44T, so that it 
is not difficult to extrapolate the value of the energy to 
infinite time. If we use the computation of the interface 



energy done by using SRSB theory |12| we expect that 
A(T) = 2.5A(T), which is very well satisfied by our data. 

While we find that the equality (^) is correct above and 
at the critical temperature (with less that a relative 1% 
error), it is definitely violated below T c ; at T = 0.7 we 
find C(l) = 0.612 ± .001 and G(l) = 0.56 ± 0.01, at T = 
0.35 we find C(l) = 0.802±0.001 and G(l) = 0.67±0.01. 
The failure of the equality ^ implies the existence of dif- 
ferent ergodic components. The q — q correlation function 
depends on the choice of the component, in agreement 
with the main prediction of the SRSB theory. 

As a last evidence we discuss the value of the local 
overlap in a box of side R, qii(x) = R~ D Yly^x+yTx+y, 
where y is an integer vector which takes all the R D val- 
ues compatible with the conditions < y v < R. We 
evaluate the probability distribution Pr^r) of the local 
box overlap. In the mean field SRSB limit the function 
Pr(qr) is Gaussian, but in a finite (not too large) num- 
ber of dimensions it is quite natural to expect deviations 
from the Gaussian limit. 

On the contrary in the droplet model the function 
Pr(ir) should have two peaks at qR ~ qEA, and should 
become the sum of two delta functions in the limit 
R — > oo (at least for R -C L). Indeed here the quan- 
tity qR is different from qEA with a probability that goes 
to zero as a power of R. 

We have focused on the Binder parameter g(R,t) = 

I — J'?!, measured after t Monte Carlo steps. At a 
given temperature we expect that for large R the data 
will collapse as 

g(R,t)=f((mr 1 ) S ) ■ (7) 

In Fig. (||) we show data for T = 0.7 from R — 3 and 
4. The scaling law we are proposing works very well. 
The Binder function extrapolate to something definitely 
different from 1, which would be the prediction of the 
droplet model, since in that case the distribution should 
be asymptotically a pair of delta function. The data ob- 
tained with the alternative temperature scheduling, by 
relatively slow cooling, give similar results for the ex- 
trapolated value of the Binder cumulant. 

We can summarize by saying that none of our find- 
ings give support to the predictions of the droplet model, 
while the broken replica approach is able to predict qual- 
itatively (and in a few cases even quantitatives) the be- 
havior of the 3c? spin glass systems. 
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FIG. 1. log (G(x, oo)) as a function of log (a;). The upper 
line from cooling, the lower one from the normal dynamics 
(see the text for details). 
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FIG. 2. The logarithm of the Binder parameter for the 
box overlap versus rescaled ratio of time and distance. Stars 
are for R = 2, hexagons are for R — 3 and asterisks for R = 4. 
The straight line is only to guide the eye. 
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